/* ****************************************** */
/* * Black Scholes Model - Parallel Version * */
/* ****************************************** */

#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>

#include "../include/rand.h"

#define MAX(a,b) ((a) > (b) ? a : b)
#define NUM_THREADS 8

/*Input variables*/
double s,e,r,sigma,T,m;

/*Trials vector*/	
double **trialVec;

/*Random numbers generator*/
struct BoxMullerState state[NUM_THREADS];

/*Size of vector of trials*/
int sizeVec;

/*Method to be executed by the threads*/
void *blackScholes(void *threadIdArg){
	int *threadIdPtr = (int *)threadIdArg;
	int threadId = *threadIdPtr;
	int i;
	double trial,aux = 0,t;

	/*Calculate trials*/
	for(i = 0; i < sizeVec; i++){
		t = s*exp((r - 0.5*pow(sigma,2))*T + sigma*sqrt(T)*boxMullerRandom(&state[threadId]));
		aux = t-e;
		trial = exp(-r*T)*MAX(aux,0);
		trialVec[threadId][i] = trial;		
	}

	/*Exit thread*/
	pthread_exit(NULL);
}

int main(){
	
	int i,j, threadError,threadId[NUM_THREADS];	
	void * status;
	
	/*Output variables*/
	double mean,stddev,confwidth,confmin,confmax;
	
	/*Threads*/
	pthread_t threads[NUM_THREADS];

	/*Standard deviation*/
	double sdAuxSum,sdAuxSumSqr;

	mean 		= 0;
	sdAuxSum 	= 0;
	sdAuxSumSqr 	= 0;
	confmin 	= 0;
	confmax		= 0;
	confwidth 	= 0;

	/*Initialize seed*/
	srand(time(NULL)*rand());

	/*Read input parameters*/
	scanf("%lf",&s);
	scanf("%lf",&e);
	scanf("%lf",&r);
	scanf("%lf",&sigma);
	scanf("%lf",&T);
	scanf("%lf",&m);

	/*Allocate result vector*/
	trialVec = (double **)malloc(NUM_THREADS*sizeof(double*));
	sizeVec  =  m/NUM_THREADS;
	for(i = 0; i < NUM_THREADS ; i++){
		trialVec[i] = (double *)malloc(sizeVec*sizeof(double));
		initBoxMullerState(&state[i]);		
	}
	
	/*Initialize threads*/
	for(i = 0; i < NUM_THREADS; i++){
		threadId[i] = i;
		threadError = pthread_create(&threads[i],NULL,blackScholes,(void *)&threadId[i]);
		
		if(threadError){
			printf("Error in pthread_create\n");
			exit(-1);
		}
	}
	
	/*Link threads to the main thread*/
	for(i = 0; i < NUM_THREADS; i++){
		threadError = pthread_join(threads[i],&status);
		if(threadError){
			printf("Error in pthread_join %ld\n",(long)status);
			exit(-1);
		}
	}
	
	/*Calculate standard deviation*/
	for(i = 0; i < NUM_THREADS; i++){
		for(j = 0; j < sizeVec; j++){
			mean += trialVec[i][j];
			sdAuxSumSqr+= pow(trialVec[i][j],2);
		}
	}
	
	/*Calculate mean, minimum and maximum confidence*/
	sdAuxSum = pow(mean,2);
	stddev = sqrt((m*sdAuxSumSqr - sdAuxSum )/(m*(m-1)));
	confwidth = (1.96*stddev)/sqrt(m);
	mean/=m;
	confmin = mean - confwidth;
	confmax = mean + confwidth;

	/*Free memory*/
	for(i = 0; i < NUM_THREADS; i++){
		free(trialVec[i]);
	}	
	free(trialVec);

	/*Output result*/
	printf("S\t%.0f\nE\t%.0f\nr\t%.0f\nsigma\t%.0f\nT\t%.0f\nM\t%.0f\n",s,e,r,sigma,T,m);
	printf("Confidence interval: (%.4lf, %.4lf)\n",confmin,confmax);

	pthread_exit(NULL);
	return 0;
}
